home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-01 / ddj9304.zip / STEREO.ZIP / LISTING1.C next >
INI File  |  1993-01-18  |  24KB  |  634 lines
  1. [LISTING1.C]
  2.  
  3. #include <file.h>
  4. #include <stdio.h>
  5. #include <unixio.h>
  6. #include <stdlib.h>
  7. #include <math.h>
  8. #include <string.h>
  9.  
  10. #define BOOLEAN     int
  11. #define TRUE          1
  12. #define FALSE         0
  13. #define OK            0
  14. #define NOT_OK       -1
  15. #define WHITE         1
  16. #define LEFT_EYE      0
  17. #define RIGHT_EYE     1
  18. #define INTEROCULAR_DISTANCE  2.5     /* inches */
  19. #define PIXEL_WIDTH           0.021   /* inches */
  20. #define ACCEPT     TRUE
  21. #define REJECT    FALSE
  22. #define SCREEN_HEIGHT         512
  23. #define NUM_DIMENSIONS          3
  24. #define WINDOW_BORDER_WIDTH     3
  25. #define NUM_SHADES            256
  26. #define MAX_LINE_LENGTH       256
  27. #define MAX_IMAGE_SIZE        200      /* 400 x 400 samples */
  28. #define MAX_LINE_SIZE         400      /* should be set to max_image_size
  29.              if there is enough main memory to hold the whole image. */
  30.  
  31. #define ROUND(x)   ((x) > 0.0  ?  (int)((x) + 0.5) : (int)((x) - 0.5))
  32. #define MAX(A, B)  ((A) > (B)  ?  (A) : (B))
  33. #define MIN(A, B)  ((A) < (B)  ?  (A) : (B))
  34.  
  35. #define LEX   TRUE
  36. #define PC   FALSE
  37.  
  38. struct    point_3D
  39. {
  40.     float  x, y, z;
  41. };
  42. typedef    struct    point_3D    point_3D_t;
  43.  
  44. struct    point_3D_ex
  45. {
  46.     float   x, y, z;
  47.     float   sha, shb;        /* shading values for triangle A and B */
  48. };
  49. typedef    struct    point_3D_ex    point_3D_ex_t;
  50.  
  51. extern int num_samples;        /* number of samples per scan line */
  52. extern int log_input;          /* 0 if data is linear, 1 if log */
  53. extern float z_scale;          /* vertical scale */
  54. extern float scan_sz;          /* 0 < X, Y < scan_sz in nm */
  55. extern point_3D_t    min, max; /* min and max of input data */
  56. extern point_3D_ex_t image[ MAX_IMAGE_SIZE ][ MAX_IMAGE_SIZE ];
  57.  
  58. /* Clipping rectangle boundaries - accessable to all routines */
  59. extern double  y_bottom, y_top, x_right, x_left, z_front, z_back;
  60. extern double  y_bottom_d, y_top_d, x_right_d, x_left_d, z_front_d, z_back_d;
  61.  
  62. extern point_3D_t  eye_pt;
  63. extern point_3D_t  light_source;
  64. extern point_3D_t  vup;
  65. extern point_3D_t  vpn;
  66. extern double      proj_plane;
  67.  
  68. /* some global variables */
  69. extern int  f_color;
  70.  
  71. /* Local variables */
  72. static double cos_theta, sin_theta;   /* to speed up rotation, dramatically! */
  73.  
  74. #if PC
  75. /*------------------------------------------------------------------------
  76.   LEX/90 compatible line drawing command.
  77. --------------------------------------------------------------------------*/
  78. void dsvec( x0, y0, x1, y1, color )
  79. int *x0, *y0, *x1, *y1, *color;
  80. {
  81.     if ( *color > 127 )
  82.         hline( 1, *x0, *y0, *x1, *y1 );
  83.     else
  84.         hline( 0, *x0, *y0, *x1, *y1 );
  85. }
  86. #endif
  87. /*------------------------------------------------------------------------
  88.   Procedure that sets the clipping rectangle limits.
  89.   The projection plane is also set to the mid-Z coordinate.
  90. --------------------------------------------------------------------------*/
  91. void set_clip_volume( x_min, x_max, y_min, y_max, z_min, z_max )
  92. double  x_min, x_max, y_min, y_max, z_min, z_max;
  93. {
  94.     y_bottom = y_min;
  95.     y_top    = y_max;
  96.     x_right  = x_max;
  97.     x_left   = x_min;
  98.     z_front  = z_min;
  99.     z_back   = z_max;
  100. }
  101. /*------------------------------------------------------------------------
  102.   Procedure that reads the STM image data.
  103. --------------------------------------------------------------------------*/
  104. static void read_stm_image_data( fp )
  105. FILE *fp;
  106. {
  107.    register i, j;
  108.    short z_data[ MAX_LINE_SIZE ];      /* 16 bit numbers */
  109.    int count, num_lines, num_lines_loaded;
  110.    double tmp, tmp1, step, hiw;   /* half image width */
  111.    long offset;
  112.  
  113.    tmp = z_scale / 16385.0;                 /* elliminates divisions */
  114.    num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
  115.    step = scan_sz / num_samples;
  116.    hiw = (double)( step * num_lines / 2.0);
  117.    offset = (long)( 2048 );
  118.    if ( fseek( fp, offset, 0 ) != 0 )
  119.       printf( "Failed to skip the STM data file header!!!\n" );
  120.    num_lines_loaded = 0;
  121.    for( i = 0; i < num_lines; i++ )   /* assumes a square image */
  122.    {
  123.       /* Read the whole line even though some of it may be thrown away. */
  124.       count = fread( (char *)( z_data ), sizeof( short ), num_samples, fp );
  125.       if ( count != num_samples )                 /* reached EOF */
  126.          break;
  127.       else  num_lines_loaded++;
  128.  
  129.       /* Assume the data is linear, for now. */
  130.       for( j = 0; j < num_lines; j++ )
  131.          image[ i ][ j ].z = (float)(z_data[ j ]) * tmp;
  132. #if 0
  133.       else if ( log_input == 1 )
  134.       {                           /* data is logarithmic (base e) */
  135.          tmp1 = 1.0 / 10.0;       /* 10 angstroms in 1 nanometer. */
  136.          for( j = 0; j < num_lines; j++ )
  137.          {
  138.             /* When STM data is logarithmic, I believe that the magnitude
  139.                is logarithmic and then the sign is added.  Thus, to expontiate
  140.                properly we must expontiate the absolute value of the data. */
  141.             if ( z_data[ j ] < 0 )
  142.                image[ i ][ j ].z = (float)( -exp( (double)(-z_data[j])
  143.                                                   * tmp ) * tmp1 );
  144.             else
  145.                image[ i ][ j ].z = (float)( exp( (double)(z_data[j])
  146.                                                  * tmp ) * tmp1 );
  147.          }
  148.       }
  149. #endif
  150.       /* Center the image around the origin, in the X and Y dimensions.
  151.          Setup X such that dx > 0, setup Y such that dy < 0. */
  152.       for( j = 0; j < num_lines; j++ )
  153.       {
  154.          image[ i ][ j ].x = (float)( j * step - hiw );
  155.          image[ i ][ j ].y = (float)( hiw - i * step );
  156.       }
  157.    }
  158.    max.x = min.x = image[ 0 ][ 0 ].x;
  159.    max.y = min.y = image[ 0 ][ 0 ].y;
  160.    max.z = min.z = image[ 0 ][ 0 ].z;
  161.    for( i = 0; i < num_lines; i++ )
  162.       for( j = 0; j < num_lines; j++ )
  163.       {
  164.          if ( image[ i ][ j ].x > max.x )   max.x = image[ i ][ j ].x;
  165.          if ( image[ i ][ j ].x < min.x )   min.x = image[ i ][ j ].x;
  166.          if ( image[ i ][ j ].y > max.y )   max.y = image[ i ][ j ].y;
  167.          if ( image[ i ][ j ].y < min.y )   min.y = image[ i ][ j ].y;
  168.          if ( image[ i ][ j ].z > max.z )   max.z = image[ i ][ j ].z;
  169.          if ( image[ i ][ j ].z < min.z )   min.z = image[ i ][ j ].z;
  170.       }
  171.    printf( "%f < X < %f\n", min.x, max.x );
  172.    printf( "%f < Y < %f\n", min.y, max.y );
  173.    printf( "%f < Z < %f\n", min.z, max.z );
  174.    printf( "Loaded %d scan lines of image data.\n", num_lines_loaded );
  175. }
  176. /* Procedure to find MIN/MAX of the image. */
  177. find_min_max( s )
  178. point_3D_ex_t  s[][ MAX_IMAGE_SIZE ];
  179. {
  180.    register  i, j, num_lines;
  181.  
  182.    num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
  183.    max.x = min.x = s[ 0 ][ 0 ].x;
  184.    max.y = min.y = s[ 0 ][ 0 ].y;
  185.    max.z = min.z = s[ 0 ][ 0 ].z;
  186.    for( i = 0; i < num_lines; i++ )
  187.       for( j = 0; j < num_lines; j++ )
  188.       {
  189.          if ( s[ i ][ j ].x > max.x )   max.x = s[ i ][ j ].x;
  190.          if ( s[ i ][ j ].x < min.x )   min.x = s[ i ][ j ].x;
  191.          if ( s[ i ][ j ].y > max.y )   max.y = s[ i ][ j ].y;
  192.          if ( s[ i ][ j ].y < min.y )   min.y = s[ i ][ j ].y;
  193.          if ( s[ i ][ j ].z > max.z )   max.z = s[ i ][ j ].z;
  194.          if ( s[ i ][ j ].z < min.z )   min.z = s[ i ][ j ].z;
  195.       }
  196. }
  197. /*------------------------------------------------------------------------
  198.   Procedure that reads the STM file header (ASCII).
  199. --------------------------------------------------------------------------*/
  200. static void read_stm_header( fp )
  201. FILE *fp;
  202. {
  203.     char line[ MAX_LINE_LENGTH + 2 ];
  204.  
  205.     while( fgets( line, MAX_LINE_LENGTH, fp ) != NULL )
  206.     {
  207.        /* At any time in the header, Control-Z => end of header */
  208.        if ( strchr( line, 0x1a ) != NULL )   break;
  209.        if ( strncmp( line, "data_type", strlen("data_type")) == 0 )
  210.        {
  211.           if ( strncmp( line, "data_type = 0", strlen("data_type = 0")) != 0 )
  212.           {
  213.              printf( "File data_type is not 0 => not STM data.\n" );
  214.              printf( "Can't handle non-STM data.\n" );
  215.              exit( 1 );
  216.           }
  217.        }
  218.        if ( strncmp( line, "num_samp", strlen("num_samp")) == 0 )
  219.        {
  220.           sscanf( line, "num_samp = %d", &num_samples );
  221.           if ( num_samples > MAX_IMAGE_SIZE )
  222.           {
  223.              printf( "Warning!!! The whole image can not be processed.\n" );
  224.              printf( "The image is %d x %d.\n", num_samples, num_samples );
  225.              printf( "The maximum allowable is %d x %d.\n", MAX_IMAGE_SIZE,
  226.                                                             MAX_IMAGE_SIZE );
  227.              printf( "Only the upper left portion of the image" );
  228.              printf( " will be processed.\n" );
  229.           }
  230.        }
  231.        if ( strncmp( line, "log_input", strlen("log_input")) == 0 )
  232.           sscanf( line, "log_input = %d", &log_input );
  233.        if ( strncmp( line, "z_scale", strlen("z_scale")) == 0 )
  234.           sscanf( line, "z_scale = %f", &z_scale );
  235.        if ( strncmp( line, "scan_sz", strlen("scan_sz")) == 0 )
  236.           sscanf( line, "scan_sz = %f", &scan_sz );
  237.     }
  238. }
  239. /*------------------------------------------------------------------------
  240.   Procedure that reads the STM image database.
  241. --------------------------------------------------------------------------*/
  242. read_stm_database()
  243. {
  244.     char fname[40];
  245.     FILE *fp;
  246.  
  247.     /* Prompt the user for the database file name */
  248.     printf( "Database file name? " );
  249.     scanf( "%s", fname );
  250.  
  251.     if (( fp = fopen( fname, "r" )) == NULL )
  252.     {
  253.        printf("Couldn't open file %s.\n", fname );
  254.        return( TRUE );
  255.     }
  256.     read_stm_header( fp );
  257.     printf( "num_samples = %d, log_input = %d, z_scale = %f, scan_sz = %f\n",
  258.              num_samples, log_input, z_scale, scan_sz );
  259.     fclose( fp );
  260.  
  261.     if (( fp = fopen( fname, "rb" )) == NULL )
  262.     {
  263.        printf("Couldn't open file %s.\n", fname );
  264.        return( TRUE );
  265.     }
  266.     read_stm_image_data( fp );
  267.     fclose( fp ); 
  268. }
  269. /*---------------------------------------------------------------------
  270.   Procedure to draw a border of the clipping window.  The border is
  271.   completely on the inside of the window.
  272. -----------------------------------------------------------------------*/
  273. void draw_clip_window()
  274. {
  275.     register  i;
  276.     int  x_min, x_max, y_min, y_max, color;
  277.  
  278.     color = 255;
  279.     x_min = ROUND( x_left );
  280.     x_max = ROUND( x_right );
  281.     y_min = ROUND( y_bottom );
  282.     y_max = ROUND( y_top );
  283.     for( i = 0; i < WINDOW_BORDER_WIDTH; i++ )
  284.     {
  285.         dsvec( &x_min, &y_min, &x_max, &y_min, &color );
  286.         dsvec( &x_max, &y_min, &x_max, &y_max, &color );
  287.         dsvec( &x_max, &y_max, &x_min, &y_max, &color );
  288.         dsvec( &x_min, &y_max, &x_min, &y_min, &color );
  289.         x_min++;  y_min++;  x_max--;  y_max--;
  290.     }
  291. }
  292. /*---------------------------------------------------------------------
  293.   Procedure that multiplies two matrices together.
  294.   Assumes that memory has been proper
  295. -----------------------------------------------------------------------*/
  296. matrix_mult( m1, r1, c1, m2, r2, c2, result )
  297. double m1[];
  298. int r1, c1;      /* number of rows and columns in matrix 1 */
  299. double m2[];
  300. int r2, c2;
  301. double result[];  /* resultant matrix */
  302. {
  303.    register i, j, k;
  304.    double c;
  305.  
  306.    if ( c1 != r2 )     /* Can't be done! */
  307.       return( 1 );
  308.    /* The resultant matrix is (r1 x c2).  So, make every element in the
  309.       resultant matrix - one element at a time. */
  310.    for( i = 0; i < r1; i++ )
  311.       for( j = 0; j < c2; j++ )
  312.       {
  313.          for( c = 0.0, k = 0; k < c1; k++ )
  314.             c += m1[ i * c1 + k ] * m2[ k * c2 + j ];
  315.          result[ i * c2 + j ] = c;
  316.       }
  317.    return( 0 );
  318. }
  319. /*------------------------------------------------------------------------
  320.   Procedure that determines the average value in the image data.
  321. --------------------------------------------------------------------------*/
  322. float average()
  323. {
  324.    register i, j;
  325.    int  num_lines;
  326.    double sum;
  327.  
  328.    sum = 0.0;
  329.    num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
  330.    for( i = 0; i < num_lines; i++ )
  331.       for( j = 0; j < num_lines; j++ )
  332.          sum += image[ i ][ j ].z;
  333.    return( (float)( sum / ( (double)num_samples * (double)num_samples )));
  334. }
  335. /*------------------------------------------------------------------------
  336.   Procedure that displays the vertex database as lines between successive
  337.   vertices.
  338. --------------------------------------------------------------------------*/
  339. display_stm_database()
  340. {
  341.     short x, y, num_lines, intensity;
  342.     float av;
  343.     short  arg[ MAX_IMAGE_SIZE * 3 ];
  344.     double z_range;
  345.     int i;
  346.  
  347.     z_range = max.z - min.z;
  348.     num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
  349.     for( y = 0; y < num_lines; y++ )
  350.     {
  351.        for( x = i = 0; x < num_lines; x++ )
  352.        {
  353.           intensity = ROUND(( image[ y ][ x ].z - min.z ) / z_range
  354.                       * (double)( NUM_SHADES ));
  355.           arg[i++] = x;
  356.           arg[i++] = y;
  357.           arg[i++] = intensity;
  358.        }
  359.        dsrnw( &num_lines, arg );
  360.     }
  361.     return( OK );
  362. }
  363. /*-------------------------------------------------------------------
  364.   Procedure that lets the user set the light source direction vector.
  365. ---------------------------------------------------------------------*/
  366. void set_light_source()
  367. {
  368.    point_3D_t ls;
  369.    double d;
  370.  
  371.    printf( "Present normalized light source vector is " );
  372.    printf( "%f, %f, %f\n", light_source.x, light_source.y, light_source.z );
  373.    printf("Specify light source direction vector from point to the origin.\n");
  374.    printf( "Enter x, y & z: " );
  375.    scanf( "%f%f%f", &(ls.x), &(ls.y), &(ls.z));
  376.    if ( ls.z < 0.0 )
  377.    {
  378.       printf( "Warning!  Z must be > 0\n" );
  379.       printf( "Light source direction unchanged.\n" );
  380.       return;
  381.    }
  382.    /* Make it a unit vector to allow easier interpretation of dot product
  383.       results. */
  384.    d = sqrt((double)( ls.x * ls.x + ls.y * ls.y + ls.z * ls.z ));
  385.    light_source.x = (float)( ls.x / d );
  386.    light_source.y = (float)( ls.y / d );
  387.    light_source.z = (float)( ls.z / d );
  388.    printf( "Negative normalized light source vector " );
  389.    printf( "%f, %f, %f\n", light_source.x, light_source.y, light_source.z );
  390. }
  391. /*-------------------------------------------------------------------
  392.   Procedure that lets the user set the viewer eye point.
  393. ---------------------------------------------------------------------*/
  394. void set_viewer_eye_pt()
  395. {
  396.    printf( "Present viewer eye point is " );
  397.    printf( " %f, %f, %f\n", eye_pt.x, eye_pt.y, eye_pt.z );
  398.    printf( "Specify viewer eye point.\n");
  399.    printf( "The viewing direction will be from this point to origin.\n" );
  400.    printf( "Enter x, y & z: " );
  401.    scanf( "%f%f%f", &(eye_pt.x), &(eye_pt.y), &(eye_pt.z));
  402.    printf( "Viewer eye point is" );
  403.    printf( " %f, %f, %f\n", eye_pt.x, eye_pt.y, eye_pt.z );
  404.    /* Set the view plane normal vector (from the eye_pt to origin). */
  405.    /* Vectors are defined as arrow point minus tail point. */
  406.    vpn.x = -(eye_pt.x);   vpn.y = -(eye_pt.y);  vpn.z = -(eye_pt.z);
  407.    printf( "View plane normal vector is" );
  408.    printf( " %f, %f, %f\n", vpn.x, vpn.y, vpn.z );
  409. }
  410. /*-------------------------------------------------------------------
  411.   Procedure that lets the user set the view up vector to other than default.
  412. ---------------------------------------------------------------------*/
  413. void set_view_up_vector()
  414. {
  415.    printf( "View plane normal vector is %f, %f, %f\n", vpn.x, vpn.y, vpn.z );
  416.    printf( "Present view up vector is " );
  417.    printf( "%f %f %f\n", vup.x, vup.y, vup.z );
  418.    printf( "Specify view UP vector from the origin to a point.\n" );
  419.    printf( "Enter x, y & z: " );
  420.    scanf( "%f%f%f", &(vup.x), &(vup.y), &(vup.z));
  421. }
  422. /*-------------------------------------------------------------------
  423.   Procedure that lets the user set the perspective projection plane
  424.   on the negative Z-axis to other than default.
  425. ---------------------------------------------------------------------*/
  426. void set_per_projection_plane()
  427. {
  428.    double e_v, e_w;     /* interocular distance in the view and world
  429.                            coordinate systems. */
  430.    double d;
  431.  
  432.    printf( "Present perspective projection plane is %lf\n", proj_plane );
  433.    printf( "For optimum viewing distance of 18 inches from the screen\n" );
  434.    e_v = INTEROCULAR_DISTANCE / PIXEL_WIDTH;
  435.    e_w = e_v / (( x_right_d - x_left_d ) / ( x_right - x_left ));
  436.    d = 18.0 * e_w / 2.5;     /* from p. 82 of Dr. Hodges' thesis */
  437.    printf( "The projection plane should be set to %lf\n", -d );
  438.    printf( "Specify new perspective projection plane on neg. Z-axis.\n" );
  439.    printf( "Enter z: " );
  440.    scanf( "%lf", &(proj_plane));
  441. }
  442. /*-------------------------------------------------------------------
  443.   Procedure that breaks the image into triangles and computes the dot
  444.   product between the normal of each triangle and the light source vector.
  445.   This information is necessary for subsequent shading computations.
  446.   The normal to a triangle defined by P11, P21, & P22 is P11P21 x P11P22.
  447.   The stepping distance in the X and Y dimensions is known and is constant,
  448.   as the image data has not been transformed yet (assumption).
  449.   Computations are minimized by removing redundant calculations and
  450.   relying on a uniform grid (dx == dy).
  451. ---------------------------------------------------------------------*/
  452. void compute_normals()
  453. {
  454.    register i, j;
  455.    point_3D_t p11_p21, p11_p22, p11_p12, cross_a, cross_b;
  456.    double d, dx;              /* stepping distance (dx = dy) */
  457.    double dx_sqrd;            /* dx^2 */
  458.    int num_lines;
  459.    point_3D_t na, nb;         /* normals to triangle A and B */
  460.  
  461.    dx = scan_sz / num_samples;
  462.  
  463.    num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
  464.    num_lines--;
  465.    dx_sqrd = dx * dx;
  466.    for( i = 0; i < num_lines; i++ )
  467.       for( j = 0; j < num_lines; j++ )
  468.       {
  469. #if DEBUG
  470.          printf( "P11 %f %f %f\n", image[i][j].x, image[i][j].y,
  471.                                    image[i][j].z );
  472.          printf( "P21 %f %f %f\n", image[i+1][j].x, image[i+1][j].y,
  473.                                    image[i+1][j].z );
  474.          printf( "P12 %f %f %f\n", image[i][j+1].x, image[i][j+1].y,
  475.                                    image[i][j+1].z );
  476.          printf( "P22 %f %f %f\n", image[i+1][j+1].x, image[i+1][j+1].y,
  477.                                    image[i+1][j+1].z );
  478. #endif
  479.          p11_p21.z = image[ i + 1 ][ j ].z - image[ i ][ j ].z;
  480.          p11_p22.z = image[ i + 1 ][ j + 1 ].z - image[ i ][ j ].z;
  481.          p11_p12.z = image[ i ][ j + 1 ].z - image[ i ][ j ].z;
  482. #if DEBUG
  483.          printf( "dz11_21 = %f,  dz11_22 = %f,  dz11_12 = %f\n",
  484.                  p11_p21.z, p11_p22.z, p11_p12.z );
  485. #endif
  486.          /* It's possible to eliminate one more multiplication in the
  487.             computations below. */
  488.          na.x = dx * ( p11_p21.z - p11_p22.z );
  489.          na.y = dx * p11_p21.z;
  490.          na.z = dx_sqrd;
  491. #if DEBUG
  492.          printf( "Na %f %f %f\n", na.x, na.y, na.z );
  493. #endif
  494.          nb.x = (-dx) * p11_p12.z;
  495.          nb.y = dx * ( p11_p22.z - p11_p12.z );
  496.          nb.z = dx_sqrd;
  497. #if DEBUG
  498.          printf( "Nb %f %f %f\n", nb.x, nb.y, nb.z );
  499. #endif
  500.          /* Normalize the normal vectors, since the intensity will be
  501.             proportional to the angle between light source and the normal. */
  502.          d = sqrt((double)( na.x * na.x + na.y * na.y + na.z * na.z ));
  503.          na.x = na.x / d;
  504.          na.y = na.y / d;
  505.          na.z = na.z / d;
  506. #if DEBUG
  507.          printf( "Na %f %f %f\n", na.x, na.y, na.z );
  508. #endif
  509.          d = sqrt((double)( nb.x * nb.x + nb.y * nb.y + nb.z * nb.z ));
  510.          nb.x = nb.x / d;
  511.          nb.y = nb.y / d;
  512.          nb.z = nb.z / d;
  513. #if DEBUG
  514.          printf( "Nb %f %f %f\n", nb.x, nb.y, nb.z );
  515. #endif
  516.          /* Compute the dot product between the light source vector and
  517.             the normals (== to the angle between two unit vectors ).
  518.             -1 <= cos( theta ) <= 1, which will be very useful. */
  519.          image[ i ][ j ].sha = light_source.x * na.x + light_source.y * na.y +
  520.                                light_source.z * na.z;
  521.          image[ i ][ j ].shb = light_source.x * nb.x + light_source.y * nb.y +
  522.                                light_source.z * nb.z;
  523.       }
  524. }
  525. /*-------------------------------------------------------------------
  526.   Sets a square matrix to identity.
  527. ---------------------------------------------------------------------*/
  528. void  set_to_identity( m, n )
  529. double m[];
  530. int n;          /* n x n matrix */
  531. {
  532.    register  i, j;
  533.  
  534.    for( i = 0; i < n; i++ )
  535.       for( j = 0; j < n; j++ )
  536.          if ( i == j )   m[ i * n + j ] = 1.0;
  537.          else            m[ i * n + j ] = 0.0;
  538. }
  539. /*-------------------------------------------------------------------
  540.   Copies the source matrix to the destination matrix of the same size.
  541. ---------------------------------------------------------------------*/
  542. void  copy_matrix( s, d, r, c )
  543. double s[], d[];
  544. int r, c;
  545. {
  546.    register i, j;
  547.  
  548.    for( i = 0; i < r; i++ )
  549.       for( j = 0; j < c; j++ )
  550.          d[ i * c + j ] = s[ i * c + j ];
  551. }
  552. /* Procedure that prints out a r x c matrix. */
  553. void show_matrix( m, r, c )
  554. double m[];
  555. int r, c;
  556. {
  557.    register i, j;
  558.  
  559.    printf( "\n" );
  560.    for( i = 0; i < r; i++ )
  561.    {
  562.       for( j = 0; j < c; j++ )
  563.          printf( "%lf\t", m[ i * c + j ] );
  564.       printf( "\n" );
  565.    }
  566. }
  567. /* Procedure to clip the ready to render polygons to the viewport.
  568.    Returns ACCEPT if the tirangle is completely  inside the viewport, or
  569.            REJECT if the triangle is completely outside the viewport. */
  570. clip_to_viewport( v, n, y )
  571. short v[];
  572. int n, y;
  573. {
  574.    register i;
  575.  
  576.    for( i = 0; i < n; i += 2 )
  577.    {
  578.       if (( v[i]   < x_left_d   ) || ( v[i]   > x_right_d ))
  579.          return( REJECT );
  580.       if ( y == LEFT_EYE )
  581.       {
  582.          if (( v[i+1] < y_bottom_d ) || ( v[i+1] > y_top_d   ))
  583.             return( REJECT );
  584.       }
  585.       else if ( y == RIGHT_EYE )
  586.       {
  587.          if (( v[i+1] < ( y_bottom_d + SCREEN_HEIGHT )) ||
  588.              ( v[i+1] > ( y_top_d + SCREEN_HEIGHT  )))
  589.             return( REJECT );
  590.       }
  591.    }
  592.    return( ACCEPT );
  593. }
  594. /* Procedure to set viewport limits. */
  595. set_viewport( x_min, x_max, y_min, y_max, z_min, z_max )
  596. double x_min, x_max, y_min, y_max, z_min, z_max;
  597. {
  598.    x_left_d = x_min;
  599.    x_right_d = x_max;
  600.    y_bottom_d = y_min;
  601.    y_top_d = y_max;
  602.    z_front_d = z_min;
  603.    z_back_d = z_max;
  604. }
  605. /* Procedure to open a communication channel to the LEX/90. */
  606. open_lex()
  607. {
  608.    short  ierr, chan;
  609.  
  610.    dsopn( &ierr, &chan );
  611.    if ( ierr > 0 )
  612.    {
  613.       printf( "An error condition occured on the LEX!\n" );
  614.       exit( 1 );
  615.    }
  616.    dsclr( &(-1) );
  617.    dscsl( &2, &0, &0 );
  618.    dscer();
  619. }
  620. /* Procedure to configure the LEX/90 as a double buffer 8 bpp system.
  621.    The LUT is initialized as well. */
  622. configure_lex()
  623. {
  624.    short  arg[10], read_array;
  625.  
  626.    /* Configure the LEX as a two buffer system - one buffer displayable
  627.       at a time (640x512), 8-plane buffers and an 8x8 LUT. */
  628.    arg[0] = 2;   arg[1] = 9;
  629.    dsesc( arg, &2, read_array, &0 );
  630.    dschan( &255, &255, &255 );
  631.    dsdisp( &0, &0, &4096 );         /* enable fill */
  632.    dsllu( &0, &0, &255, &255 );     /* ramp the LUT (256 shades of gray) */
  633. }
  634.